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OQ We use the boundary element method to study the low- Reynolds number locomotion of a 



a 
>. 



spherical model microorganism in a circular tube. The swimmer propels itself by tangen- 
tial or normal surface motion in a tube whose radius is on the order of the swimmer size. 



^ Hydrodynamic interactions with the tube walls significantly affect the average swimming 



^ speed and power consumption of the model microorganism. In the case of swimming 

parallel to the tube axis, the locomotion speed is always reduced (resp. increased) for 

_ ^ swimmers with tangential (resp. normal) deformation. In all cases, the rate of work nec- 

1^,^ essary for swimming is increased by confinement. Swimmers with no force-dipoles in the 

1-^ far field generally follow helical trajectories, solely induced by hydrodynamic interactions 

'— ^ with the tube walls, and in qualitative agreement with recent experimental observations 

T-H for Paramecium. Swimmers of the puller type always display stable locomotion at a lo- 

I cation which depends on the strength of their force dipoles: swimmers with weak dipoles 

J^ (small a) swim in the centre of the tube while those with strong dipoles (large a) swim 

t — near the walls. In contrast, pusher swimmers and those employing normal deformation 

\l are unstable and end up crashing into the walls of the tube. Similar dynamics is observed 

^T") for swimming into a curved tube. These results could be relevant for the future design 

. . of artificial microswimmers in confined geometries. 
> 

K^ Key words: low-Reynolds number swimming, boundary element method, hydrodynamic 

2 interaction, swimming microorganisms 



1. Introduction 

The locomotion of self-propelled microorganisms have recently attracted sizable at- 



tention in both the applied mathematics and biophysics communities ( Lighthill 1975 
f Email address for correspondence: lailai@mech.kth.se 
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[19761 [BreSieii & Winet|T977||Purcell][T977l|Yates|T986||Berg|2000l [Fauci fc Dillon|2006[ 



Lauga fc Powers||2009). A number of novel phenomena have been discovered, including 



the dancing behaviour of pair Volvox algae (Drescher et aL|[2009l, the collective motion 



of motile Bacillus suhtilis bacteria (Dombrowski et al. 20041, and tumbling dynamics 



of flagellated Chlamydomonas ( Polin et aL||2009" Stocker fc Durhain||2009 ) . One area of 
particularly active research addresses the variation in cell mobility as a response to com- 
plex environments, including the dependence on the rheological properties of the medium 
where cells swim ( |Lauga||2007l |Fu et aL|[2008l [Elfring et ar]|2010[ |Liu et aL|[20Tl| |Shen 
fc Arratia|2011 Zhu et aL|2011 2012), the presence of an external shear flow (Hill et al. 



2007 Kaya fc Koser|20"T2), gravity (Durham et a/.|2009), or a sudden aggression (Hamel 



eFaIl[20lT| ) 



Many microorganisms swim close to boundaries, and as a result the effect of boundaries 
on fluid-based locomotion has been extensively studied. E. coli bacteria display circular 



trajectories near boundaries, clockwise when the wall is rigid (Lauga et al. 2006) and 



anti-clockwise near a free surface ( Leonardo et aL||2011 |. Experiments, simulations, and 
theoretical analysis are employed to investigate locomotion near a plane wall ( Katz|1974 
[T9751 IRamia et aLl[T993| [Fauci fc Mcdonald] [T9951 [Goto et a/.|[2005l [Berke et aL|[2008l 
Smith et a/. [2009 Shum et a/.|2010 Spagnolie fc Lauga|2012 ) explaining in particular the 



accumulation of cells by boundaries ( Ramia et a/.||1993 Fauci fc Mcdonald||1995 Berke 
et a/.|2008[ [Smith et a/.|2009[ |Shum et a/.|2010| [Drescher et al.\2QU\ . Most of these past 



studies consider the role of hydrodynamic interaction in the kinematics and energetics of 
micro-scale locomotion, developing fundamental understanding of how microorganisms 
swim in confined geometries. 

Although most past studies consider interactions with a single planar, infinite surface, 
microorganisms in nature are faced with more complex geometries. For example, mam- 



malian spermatozoa are required to swim through narrow channel-like passages (Winet 



1973 Katz 1974 ), Trypanosoma protozoa move in narrow blood vessels (Winet 1973 ), and 



bacteria often have to navigate microporous environments such as soil-covered beaches 
and river-bed sediments ( Biondi et aL|[1998 ). 

Locomotion of microorganisms in strongly confined geometries is therefore biologically 
relevant, and a few studies have been devoted to its study. An experimental investigation 



was conducted by Winet ( 1973 ) to measure the wall drag on ciliates freely swimming in a 



tube. Perturbation theory was employed to analyse the swimming speed and efficiency of 
an infinitely long model cell swimming along the axis of a tube ( Felderhof|2010 ). Numer- 
ical simulations using multiple-particle collision dynamics were carried out to study the 
motion of model microswimmcrs in a cylindrical Poiseuille flow ( Zottl fc Stark|20r2 ) . Re- 
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cent experiments ( Jana et a/.||2012 ), which originahy inspired the present paper, showed 
that Paramecium cells tend to follow helical trajectories when self-propelling inside a 
capillary tube. 

In this article, we model the locomotion of ciliated microorganisms inside a capillary 
tube. Specifically, we develop a boundary element method (BEM) implementation of the 
locomotion of the squirmer model ( Lighthill|1952 Blake|197T ) inside straight and curved 
capillary tubes. The boundary element method has been successfully used in the past 



to simulate self-propelled cell locomotion at low Reynolds numbers (Ramia et al. 1993 



Ishikawa et aL|[2006 Shum e^"fll]|2010 Nguyen et flL||20lT ). Our specific computational 
approach is tuned to deal with strong geometrical confinement whereas traditional BEM 
show inaccuracy when the tube becomes too narrow ( Pozrikidis||2005 ) . 

After introducing the mathematical model, its computational implementation and val- 
idation, we calculate the swimming speed and power consumption of spherical squirmers 
with different swimming gaits inside a straight or curved capillary tube. The effect of 
tube confinement, swimming gait, and cell position is investigated. By studying tra- 
jectories of squirmers with varying initial cell positions and orientations, we show that 
cells end up either swimming parallel to the tube axis or performing wavelike motions 
with increasing/decreasing wave magnitudes. The dynamic stability of the cell motion is 
also analysed revealing the importance of the swimming gaits. In particular, squirmers 
employing the gait leading to minimum work against the surrounding fluid are seen to 
generically execute helical trajectories, in agreement with the experimental observation 
of swimming Paramecia inside a capillary tube (Jana et aL||2012). 



2. Mathematical Model 

2.1. Squirmer model 

In this work we use steady squirming as a model for the locomotion of ciliated cells such 
as Paramecium - more specifically, as a model for the envelope of the deforming cilia 
tips at the surface of the cells. This steady model has been employed in the past to ad- 
dress fundamental processes in the physics of swimming microorganisms, such as nutrient 



uptake ( Magar et al. 2003 1 , locomotion in stratified and viscoelastic fluids ( Doostmoham 



madi et al. 2012 Zhu et al. 2012), biomixing (Lin et al. 2011 1, and the collective behaviour 



of microorganisms (Ishikawa & Pedley 2008 Underbill et al. 2008 Evans et al. 2011 ). Fur- 



thermore, simulations of two interacting Paramecium using the squirmer model showed 
good agreement with corresponding experiments ( Ishikawa fc Hota||2006 ). 

In the model, a non-zero velocity, U52-, is imposed at the surface of the spherical 
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swimmer as first proposed by Lighthill ( 1952 ) and Blake ( |1971 ). In this work, we consider 

for the most part pure tangential surface deformation (normal surface deformation will 



be covered in ^4.6 only) and adopt the concise formulation introduced in Ishikawa & 



Pedley (20081 where the imposed velocity on the surface of a squirmer centred at the 



origin is explicitly given as 



usT(r) = Y^ 



B P' 



e • r 

r 



err 
e 

r r 



(2.1) 



where e is the orientation vector of the squirmer, _B„ is the nth mode of the tangential 



surface squirming velocity (Blake 1971), P„ and P!^ are the nth Legendre polynomial 
and its derivative with respect to the argument, r is the position vector, and r = |r|. In 
a Newtonian fluid, the swimming speed of the squirmer in free space is Ugj, = 2Bi/3 
(Blake|1971) and thus dictated by the first mode only. The second mode, B2, governs the 



signature of the flow field in the far field (stresslet) . As in many previous studies ( Ishikawa 
et al. 2006 Ishikawa fc Pedley|2008 ), we assume _B„ = for n > 2. In that case, the power 
consumption by the swimmer is Vgrp = Sttim (2Bf + _B|) /3, where ^ is the dynamic 
viscosity of the fluid and a the radius of the sphere. 

The tangential velocity on the sphere in the co-moving frame is therefore simply 
expressed, in spherical coordinates, as ug{d) — Bisinfl + (i?2/2) sin26', where 6 = 
arccos(e • r/r) is the polar angle between the position vector r and the swimming di- 
rection e. We introduce an additional dimensionless parameter, a, representing the ratio 
of the second to the first squirming mode, a = B2/B1. When a is positive, the swimmer 
is called a puller and obtain the impetus from its front part. As a is negative, the cell is 
called a pusher and thrust is generated from the rear of the body. A puller (resp. pusher) 



generates jet-like flow away from (resp. towards) its sides, as shown in Ishikawa (2009) 
and references therein. A squirmer with a = is termed a neutral squirmer, and it is 
associated with a potential velocity field. 

We note that the model we consider does not capture the unsteadiness of the flow aris- 
ing from the periodic beating of flagella and cilia in microorganisms such as Paramecium 



or Volvox (Guasto et al. 2010 Drescher et al. 2010). Here we assume that the steady, 



time-averaged, velocity dominates the overall dynamics, and will consider the underlying 
unsteadiness in future work. 

2.2. Swimming in a tube 

The spherical squirmer (radius, a) is swimming in a cylindrical tube of radius R, as 
illustrated in figure [l] The centre of the squirmer is located at a distance b from the 
tube axis. We use Cartesian coordinates with an origin at the centre of the tube and 
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Figure 1 . Schematic representation of a spherical squirmer of radius a swimming in a tube of 
radius R. The centre of the squirmer is located at a distance h from the tube axis. The origin 
of the Cartesian coordinates coincides with the centre of the tube. The bounding surfaces to 
the fluid are denoted S (surface of squirmer), B (bottom tube cap), T (top tube ca,p) and C 
(surface of the tube conduit). 



the x-direction along the tube axis. As in Higdon & Muldowney ( 1995 ) we introduce the 
nondimensional position /3 as 

fi^bl{R-a), (2.2) 

so that /3 = indicates that the squirmer is at the centre of the tube while for /3 = 1 the 
squirmer is in perfect contact with the tube wall. 

3. Numerical method 

3.1. Formulation 

The boundary element method (BEM) has already been successfully adopted to study the 
hydrodynamics of swimming microorganisms in the Stokesian regime (Ramia et aL|1993 



Ishikawa et al. 2006 Shum et al. 20101. Our current work mainly follows the approach 



Pozrikidis (2002), the important difference being that we use quadrilateral elements 



instead of triangle elements as typically used and originally proposed. The method is 
introduced briefly here. 

In the Stokesian realm, fluid motion is governed by the Stokes equation 



Vp + /xV^u = 0, 



(3.1) 
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where p is the dynamic pressure and u the fluid velocity. Due to the hnearity of the Stokes 
equation, the velocity field, u(x), resulting from moving bodies with smooth boundary 
S can be expressed as 



u(x) = -^ /f(x')-S(x,x')d^x', 



(3.2) 



where f (x') is the unknown force per unit area exerted by the body onto the fluid. The 
tensor S is the Stokeslet Green's function 

S,(-x')^(| + ^). (3.3) 

with di — Xi — x'^, d^ = |x — x'p = df + d"^ + d'^, and 6ij denoting the Kronecker delta 
tensor. 

We discretize the two bodies in the problem, namely the spherical squirmer and the 
surrounding tube, into N zero-order elements with centres at the locations {x^, g = 1 — >■ 
A''}, with (7 = 1 — > Ns denoting the elements on the squirmer surface and q — Ns + l — > N 
the elements on the surface of the tube. For the r*'* element, f (x') is assumed to be 
constant over the element and is thus approximated by the value f^. As a consequence, 



the discretized version of 3.2 is, when evaluated on one of the elements, 

1 ^ r 

— Vf,-/ S(x„x')d5x',g = l^iV. (3.4) 



Stt/i 



In its discrete form, equation 3.4 represents a total of 3iV equations for the 3iV unknown 
force density components. 

3.2. Swimming and squirmer boundary conditions 

On the squirmer surface, the left-hand side of |3.4| is not fully known. The swimmer 
has an instantaneous surface deformation, ug, plus 6 unknown components, namely its 
instantaneous translational velocity vector, U, and its instantaneous rotational velocity 
vector, rj. Thus, the left hand side of|3.4| when evaluated on the surface of the squirmer, 



becomes u (x^) = U + Jl x x^ + U5 (xg) for q from 1 to Ng (here x^ = x — x-'^, where x^ 
is an arbitrary reference point, the centre of the spherical squirmer for convenience). The 
6 additional equations necessary to close the linear system are the force- and torque-free 
swimming conditions, namely 

/ f (x) dS-x = 0, / S X f (x) d5x = 0, (3.5) 

for X G squirmer. 
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3.3. Other boundary conditions 
The situation addressed in our paper is that of a squirmer swimming inside an infinitely 
long tube filled with a quiescent fluid. Numerically, we close both ends of the tube with 
appropriate boundary conditions. If the tube caps are sufficiently far away from the 
squirmer, the velocity near the caps is almost zero, so we have u^ = ut = 0, and the 
pressure over the bottom and top cap is pB and pr respectively ( Pozrikidis]|2005 1 . The 
force density fr over the top cap can be approximated by fr = prn (Pozrikidis 2005), 



where n is the unit normal vector pointing from the top cap into the fluid domain. Since 
pressure is defined up to an arbitrary constant, without loss of generality, we set px = 0, 



fr = and the top cap does not requires discretization. However, unlike Pozrikidis ( 2005 ), 



we do perform discretization on the bottom cap, solving for the normal and tangential 
components of the force density fs there. For the conduit part of the tube, we use no-slip 
boundary condition, thus write Uc — 0. 

Since we set the velocity on both caps of the tube to be zero, the error due to domain 
truncation need to be carefully considered. A truncated tube length L of 7ri? or I-kR 



was chosen in Pozrikidis (2005) and L = 3i? in Higdon & Muldowney (1995). In our 



computation of hydrodynamic force on a moving sphere inside, we tested different values 
L and examined the truncation error. We find the length, L = 2'kR, to be long enough 
for required accuracy (see figurepland details below). In the case of swimming squirmers, 
we set L ~ SttR, and larger values of L were shown to have negligible differences in the 
results. 



3.4. Discretization and integration 
Zero-order constant quadrilateral elements are used to discretize all the surfaces. We use 



six-patch structured grid to discretize the sphere (Higdon & Muldowney 1995 Cortez 



et aL|[2005 Smith|[2009 1 , mapping six faces of a cube onto the surfaces of a sphere with 
each face latticed into a square mesh. The conduit part of the tube is divided into cylindri- 
cal quadrilateral elements obtained from the intersections of evenly spaced planes normal 
to tube axis and evenly spaced azimuthal planes (Higdon fc Muldowney|19"95 Pozrikidis 



2005 Wen et al. 2007). Moreover, orange- like quadrilateral elements are used for the 



bottom cap of the tube ( Higdon fc Muldowney|1995 Wen et aL|2007 |. For the sphere we 
adopt the six-patch quadrilateral grid with parameterized coordinates instead of triangle 



elements (Pozrikidis 2002 2005). Such discretization with its natural parametrisation 



facilitates Gauss-Legendre quadrature when performing numerical integration. Template 
points used in the quadrature lie exactly on the sphere surface since their coordinates 
are derived from the parametrisation. The resulting improved quadrature gives superior 
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Figure 2. Local mesh refinement of the cyhnder (red) and the sphere (green). The geometrical 
parameters are a/R = 0.3 and (3 = 0.95. For a better visualisation, the mesh on the squirmer 
surface is reproduced on the displaced sphere as indicated by the black arrow. 

accuracy (see Table Rj). The integration for singular elements are performed by using 
plane polar coordinates with Gauss-Legendre quadrature ( Pozrikidis||2002 ) . 

In many instances, the squirmer is so close to the cylindrical wall that near-singular 
integration has to be performed, a key point to achieve the required accuracy and ef- 



ficiency (Huang & Cruse 1993). We perform local mesh refinement in the near-contact 



regions between the squirmer and the tube ( Ishikawa et aL|[2006 Ishikawa fc Hota|2006 1 
as illustrated in figure [2] The agreement between numerical results with our method 



and existing results from high-order spectral boundary element method (Higdon & Mul- 



downey|1995 ) improves significantly when applying such local mesh refinement as shown 
in the next section where we compute the resistance of a translating sphere inside a 
cylindrical tube. 

3.5. Validation and accuracy 

We first compute the drag force, F, on a translating sphere in an unbounded domain and 
compare it with the analytical expression, F = GirfiaU, where /i is the dynamic viscosity 
of the fluid and U is the translational speed of the sphere. As shown in Table [l] the cur- 
rent method is very accurate when compared to the three similar approaches ([Pozrikidis 



2002 Cortez et al. 2005 Smith 20091. We then compute the drag force and torque on 



a sphere translating parallel to an infinite, flat, no-slip surface. The surface is modelled 
by a discretized plate of size 40a x 40a. Our simulation agree well with analytical re- 
sults ( Goldman et al. 1967 ) , as shown in Table [2] Finally, we compute the drag force 
acting on a sphere translating inside the tube with confinement a/R — 0.4, up to a 
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Cortez et al. 



(20051 



Smith 



(20091 



Pozrikidis 



(20021 



Element Order 
(functional variation) 
Element Type 
Element Number 
Singular 
integration 



Quad 
6 X 12 X 12 
Regularization 
e = 0.01 



Quad 
6x6x6 
Regularization 

6 = 0.01 

with adaptive 
Gauss 
Quadrature 
0.431 

Table 1. Relative error, in percentage, in the drag force on a translating sphere in an unbounded 
domain between the method in this paper and three other methods. The parameter e is the 



Relative error (%) 



12.6 



Tri 

512 

Analytical 

integration 

with 

Gauss 

Quadrature 

9.610"^ 



This paper 




Quad 

6x6x6 

Analytical 

integration 

with 

Gauss 

Quadrature 

1.4 lO"'^ 



regularization parameter first introduced in Cortez et al. (20051. 



Ferr (%) Terr (%) 

0.00426 0.09488 

0.01911 0.37879 

0.04274 0.20478 

0.07809 0.25773 

0.09405 0.74217 

1.005004 0.17669 1.13493 

1.003202 0.27472 1.74313 

Table 2. Relative error in the drag force, F^rr, and torque, T^rr, in percentage, on a sphere 



h/a 
3.7622 
2.3523 
1.5431 
1.1276 
1.0453 



translating parallel to an infinite wall between our computations and the analytical results ( Gold- 



man et ar||1967 l. Here a is the radius of the sphere and h the distance between the centre of 
the sphere and the wall. 



maximuni value oi f3 — 0.99, and compare our results with published data obtained with 



high-order spectral boundary element method (Higdon & Muldowney 19951. As illus 



trated in figure |3J the maximum relative error is less than 1.2%. In all simulations, the 
maximum confinement is taken to be /3 = 0.99 to ensure sufficient accuracy. 



10 



Rp 



Lailai Zhu, Eric Lauga and Luca Brandt 
(%) 




0.10.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



Figure 3. Relative error in the three components of the drag force, R^^, Ry, and -R^, on a sphere 



translating inside a tube {a/R = 0.4) between the present paper and Higdon & Muldowney 



( 1995 I. Note that the three largest values of /3 chosen are 0.95, 0.975 and 0.99. 



4. Swimming inside a tube: results 

We now have the tools necessary to characterise the locomotion of squirmers inside 
a tube. Our computational results, presented in this section, are organised as foUovirs. 
We first compute the svirimming kinematics and power consumption of a squirmer in- 
stantaneously located at various positions inside the tube while its orientation is kept 
parallel to the tube axis. These results then enable us to understand the origin of the 
two-dimensional wave-like trajectory for a neutral squirmer inside the tube. We also anal- 



yse the asymptotic stability of trajectories close to solid walls (Yizhar fc Richard||2009 



Crowdy & Yizhar 2010). We then move on to examine the general three-dimensional 



helical trajectory of a neutral squirmer and also consider the kinematics of pusher and 
puller swimmers. Finally, we study locomotion induced by normal surface deformation 
and consider locomotion inside a curved tube. 

4.1. Static kinematics and energetics 

To start our investigation, we first numerically calculate the swimming speed and power 
consumption for a squirmer exploiting pure tangential surface deformation (for complete- 



ness, results on squirmers with normal surface deformations are shown in Sec. 4.6 1. We fix 



Bi — 1 and vary the value of a, while different values of a/R and f3 are chosen to address 
the effect of confinement and eccentricity on the instantaneous swimming kinematics. 

In figure |4J we plot the instantaneous swimming speed of a squirmer with orientation 
parallel to the tube axis (positive x direction) and location (a;, y, z) — (0, 0, — /3 {R — a)). 
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Figure 4. Instantaneous swimming speed for a squirmer in a capillary tube. Left: Swimming 
velocity in the axial direction, [/x (independent of a), scaled by the swimming speed in free space 
Ust'^ different values of a/R are reported with a maximum value of /3 = 0.99. Right: Swimming 
velocity Uz in the transverse direction, scaled by the swimming speed in free space; here a/R is 
fixed to 0.3. Different values of a are reported with maximum value of /3 = 0.99. In both figures, 
the squirmer is located at (0,0, — /3 {R — a)) with its orientation parallel to the tube axis. 



The swimming velocity parallel to the tube axis {Ux) is displayed in figure [4] (left) while 
the velocity perpendicular to it (Uz) is shown in figure l4] (right). Interestingly, both 
pushers and pullers have the same swimming speed, Ux, as the neutral squirmer. This 
is due to the fact that the second squirming mode, ^ B2 sin 29, is front back symmetric, 



and thus produces zero wall-induced velocity (Berke et al. 2008), as confirmed by our 



simulation. We observe numerically that when a = 0, there is only one non-zero velocity 
component, namely Ux- In contrast, for pushers and pullers (a 7^ 0) a non-zero transverse 
velocity component, Uz, is induced. The value of Ux is seen to decrease with confinement, 
a/R, and eccentricity, /3, as shown in figured (left). The sharp decrease when /3 is beyond 
Ri 0.8 is due to the strong drag force experienced closer to the wall which overcomes the 
propulsive advantage from near-wall locomotion. 

The transverse velocity, Uz, shown in figure El (right), is plotted against the swimmer 
position, /3, for different values of a while the confinement is fixed at a/R = 0.3. In the 
case of a puller [a > 0), the swimmer will move away from the nearest wall [Uz > 0) 
while a pusher (a < 0) will move towards the nearest boundary [Uz < 0), as expected 



considering the dipolar velocity field generated by squirmers (see also Sec. 4.4 1. The abso- 



lute value of of Uz increases with a and is of the same magnitude for pushers and pullers 



of equal and opposite strength. A similar effect was explained in Berke et al. ( 2008 1 for 



a plane wall, although in that case, the cell was approximated by a point stresslet and 
the cell-wall distance was considerably larger than the cell size. By probing hydrodynam- 
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Figure 5. Left: Rotational velocity of the squirmer in the direction normal to the plane of 
locomotion, Qy. The squirmer is located at (0, 0, —(3 (R — a)) with its orientation parallel to the 
tube axis. Different values of a/R are plotted with maximum value oi P — 0.99. Right: Physical 
picture of cell rotation near the walls. The circle indicates the spherical squirmer and the green 
arrows denote the tangential surface velocity imposed by the squirmer. Blue straight arrow and 
curved arrow denote the cell orientation and rotational velocity respectively; 1 indicates the 
closest point on the cell to the top wall and 2 the closest point to the bottom wall while F/" 
and Fi' are the wall-induced shear forces generated near point 1 and 2. 



ics very close to the wall, we observe that the magnitude of Uz does actually not vary 
monotonically with /3, instead reaching a maximum value as /? « 0.9. Moving away from 
the tube centre, the transverse velocity increases due to stronger hydrodynamic interac- 
tions with the tube walls before decreasing owing to a significantly larger hydrodynamic 
resistance very close to the tube boundaries. 

Beyond the translational velocities, the squirmers also rotate due to hydrodynamic 
interactions with the tube boundaries. Numerical results show that the magnitude of the 
rotational velocity, Jl, is independent on the dipole strength, a, and that all squirmers 
rotate away from the closest wall. This is also attributed to the front-back symmetric 
distribution of the second squirming mode. Using our notation, we therefore obtain that 
squirmers rotate in the ~y direction. The value of fiy is displayed in figure p^ (left). Its 
magnitude increases with eccentricity, /3, and confinement, a/R, as a result of stronger 
hydrodynamic interactions. To explain the sign of the rotational velocity, we look in 
detail at a neutral squirmer in figure [S] (right), in the case where the swimmer is located 
closer to the top wall. Green arrows display the tangential surface deformation which 
generates locomotion. Given points 1 and 2 on the squirmer surface, the black arrows 
indicate the velocity field and show that the shear rate is higher near point 1 than point 
2. Consequently, the wall-induced force on point 1, Ff, is larger than that on point 2, 
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Figure 6. Power consumption, V , of the neutral squirmer (a = 0, left) and puller (a = 5, right), 
scaled by their corresponding values in free space. Insets display P as a function of ajR in the 
case /3 = 0. The orientation of the squirmer is parallel to the tube axis with the maximum value 
of /3 = 0.99. 



F2 1 producing a resultant clockwise torque. Since the total torque on the squirmer is 
zero, the squirmer has to rotate in the clockwise direction to balance this torque, escaping 
from the top wall. When the squirmer is closer to the top wall, an increased asymmetry 
will induce a stronger rotation. 

Next, we analyse the power consumption by the squirmer. The power, V ^ is defined as 
'P = Js iout ■ ugrfS', where tout is the force per unit area exerted from the outer surface 
of the body onto the fluid and U5 is the squirming velocity. In the single-layer potential 

as in 



formulation, equation 3.2 



Ishikawa et al. ( 2006 1 , the unknown f is the sum of the 



force density from outer (fowt) and inner (fi„) surface. We therefore rewrite the power as 
V — Jgi ■ usdS — Jg fin • usdS, where L fin ■ usdS denotes the viscous dissipation of the 
flow inside the squirmer. We thus need to subtract the internal viscous dissipation in the 
fluid given by the numerics where fin can be derived analytically based on the squirming 
velocity. In figure [6j we depict the dependence of V, scaled by the corresponding value 
in free space, with /3 and for different values of a. For each gait, V increases slowly until 
/3 Rs 0.8 followed by a rapid increase for cells closer to the wall. Such a drastic power 
increase is in agreement with the sharp decrease in swimming speed close to the tube, and 
consequently, a significant decrease in swimming efficiency is expected. In addition, as the 
confinement is getting stronger, the eccentricity of the swimmer's position becomes more 
important. For example, as /3 changes from to 0.99, the power consumption of a neutral 
squirmer V increases only by around 45% for a/R = 0.2 but by 85% for a/R = 0.5. 
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Figure 7. Sketch of the spherical squirmer inside the tube with coordinate system and angles. 
The initial dimensionless off-axis distance is measured by /?/ while the angles ^/ and xi control 
the initial cell orientation, e. For xi = 0, and in the absence of noise, the squirmer motion is 
restricted to the x — z plane. 



4.2. Two-dimensional wavelike motion of the neutral squirmer 

We next study in detail the trajectory of a squirmer inside a tube with fixed confinement; 
unless otherwise stated, all results in this section use the same value, a/R — 0.3. The cell 
is neither a pusher nor a puller, but a neutral squirmer generating potential flow field 
(a = 0). The initial position and orientation of the cell are defined as in figure l7J The cell 
is initially placed at (0,0,-6/), with bj — f3j (R— a), and oriented parallel to the axis 
(^7 = 0); the motion of the cell will also be restricted to the x ~ z plane {xi ~ 0). We 
calculate the translational and rotational velocity of the cell at each time step and update 



its position using fourth-order Adams-Bashforth scheme as in Giacche & Ishikawa (2010). 



Note that in the simulations the cell always remains in the centre of the computational 
domain (while its axial velocity is stored for post-processing), which allows to minimise 
the error introduced by domain truncation. 

Our computations show that the squirmer always displays a periodic wavelike trajec- 
tory in the tube, with amplitude A and wavelength A. This is illustrated in figure [8] for 
Pi — 0.9 (top) and f3j — 0.7 (bottom). The wave amplitude does not change over time 
and is two times the initial off-axis distance, namely, A — 26/ . The presence of a nonzero 
rotational velocity, fly, discussed above and shown in figure [51 is the key parameter lead- 
ing to the periodic trajectory. By considering cases where the initial orientation of the 
cell is not parallel to the axis (thus for which the orientation vector has non-zero x and z 
components) and we find that as long as the squirmer does not immediately descend into 
the wall, a wavelike trajectory is also obtained. To present all results in a concise man- 
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Figure 8. Two-dimensional trajectories of a neutral squirmer inside a capillary tube with con- 
finement a/_R=0.3. All positions are measured in units of cell radius a (same for figures here- 
inafter unless otherwise specified). Blue circles and red arrows indicate, respectively, the instan- 
taneous position and orientation of the squirmer. The cell is released from (0, 0, — /3/ {R — a)) 
with Pi — 0.9 (up) and /3/ — 0.7 (bottom), while the initial orientation is parallel to the axis. 
We denote A the wavelength of the periodic trajectory and A its amplitude. 



ner, we consider the motion of the neutral squirmer as a dynamical system similarly to 
recent work on two-dimensional swimming (Yizhar fc Richard||2009 Crowdy & Samson 



2011 ). The trajectory is defined by two parameters, the off-axis distance (z) and the angle 



between the swimmer orientation and the tube axis (^). We report the phase portrait 
of the neutral squirmer in the {z, ^) plane in figure ^ where the solid curves show the 
trajectories. The marginally stable point (0, 0) corresponds to locomotion along the axis 
of the tube. For any initial conditions (2,^), the neutral squirmer swims along wavelike 
trajectories corresponding to the periodic orbits in figure^ (the largest periodic orbit in 
the figure has a maximum /3 of 0.95). 

The main characteristics of the squirmers' trajectories are shown in figure[lO]for differ- 
ent initial positions, /3/. We display the trajectory wavelength, A, and the wavelength-to- 
amplitude ratio, A* = X/{A/2). It is clear that A and A* both decrease with /3/. Indeed, 
when the swimmer is at the crest or trough of the periodic trajectory, stronger rotation 
occurs for larger /?/. Therefore, the swimmer will escape from the nearest wall more 
rapidly, resulting in a decrease of the wavelength. We also show in figure [TO] that the 
time-averaged axial speed, Ux, and the time-averaged swimming speed along the trajec- 
tory, U, decrease with /?/ whereas the time-averaged power consumption, V, increases 
when the squirmers move closer to the wall. 
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Figure 9. Phase portrait for the neutral squirmer inside the tube in the (z,^) plane with 
confinement a/R = 0.3. Closed orbits correspond to 2D wavelike trajectories. The black cross 
denotes the equilibrium point {z,(^) = (0,0). 



4.3. Three-dimensional helical trajectory of the neutral squirmer 

By tilting the initial cell orientation, e, off the x ^ z plane, the squirmer trajectories 
become three dimensional and take the shape of a helix, a feature we address in this 
section. As in the two-dimensional case, these three-dimensional trajectories are a con- 



sequence of hydrodynamics interactions only. Recent experiments in Jana et al. (2012) 



showed that Paramecium cells display helical trajectories when swimming inside capillary 
tubes, a feature our simulations are thus able to reproduce. Note that some Paramecium 
cells also follow helical trajectories in free space due to asymmetries in the shape of their 
body and the beating of its cilia. In the current work we focus on swimmers whose helical 
dynamics arises only in confinement. 

We introduce xi ^^ the yaw angle between the initial cell orientation and the x — z 
plane (see figure ItI), so that the initial orientation becomes (cos (x/) , sin {xi) ,0). In our 
simulations, /3/ ranges from 0.3 to 0.9 and xi from 20° to 40°. Within these parameters, 



squirmers always display helical trajectories. One such helix is plotted in figure 11 for an 
initial position /3/ — 0.8 and a yaw angle xi — 40°. The helical trajectory is a combination 
of wavelike motions developed in the azimuthal y — z plane and in the axial direction, see 
figure [TTj j and c. In figure fTTb, we show the projected circular trajectory of the swimmer 
in the y—z plane. In figure [lip , we show that the curves y (x) and z {x) share the same 
wavelength and time period. We then plot the values of z and h (cell off-axis distance) 
as a function of the axial position, x, during one period in figure [TTJi to show that the 
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Figure 10. Dynamics and kinematics of the neutral squirmer in a tube as a function of the 
initial dimensionless off-axis position, /?/. Top left: wavelength, A, of the periodic trajectory. 
Top right: wavelength-to-amplitude ratio. A* = X/{A/2). Bottom left: time-averaged swimming 
speed in the axial direction, Ux, and along the trajectory, U, both rescaled by the free-space 
swimming velocity. Bottom right: time-averaged power consumption V, rescaled by value in free 
space. 

wave frequency of h(x) is three times that of z{x). Indeed, the trajectory projected 
in the plane perpendicular to the tube axis resembles a regular triangle (A1A2A3), 
with vertices corresponding to locations of maximum off-axis distance where the cell 
bounces back inside the tube. In this particular case, the cell bounces off the wall three 
times during one orbit with an angle ^ = 60°. A variety of other wave patterns can 
be observed for different initial cell positions and yaw angels (/3/,x/)- We display two 



of them in figure 12 in the y — z plane, with {f3i,xi) = (60,30°) (figure 12 left) and 



(60, 20°) (figure 12 right). The swimmer on the left approaches the wall 21 times during 
one periodic orbit with V' — 42.86°, whereas the example on the right displays a 5 fold 
helix with ip = 36° . 



Finally in figure 13 we show the variation of the averaged swimming speed (left) and 



power consumption (right) with the initial cell position (/?/) and orientation (xi), where 
both the speed and power are nondimensionalized by their corresponding values in free 
space. The time-averaged swimming speed along the axial direction, Ux, and along the 
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Figure 11. Three-dimensional trajectory of the neutral squirmer in the tube. The initial position 
and orientation are given by /?/ = 0.8, xi ~ 30°, ^/ = 0°. (a): trajectory in perspective view; the 
empty circle and solid triangle indicate the start and end of one periodic orbit; (b): trajectory in 
the y — z plane (axes shown in figure [lip .) ; (c): trajectories in the x — y (dashed) and x~ z (dot 
dashed) planes; (d): relation between the wavelike motion developed in the axial and azimuthal 
direction. 
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Figure 12. Periodic orbits of the neutral squirmer in the transverse plane for two different 
initial positions and orientations. Left: {Pi,Xi) ~ (0.6,30°). Right: {Pi,xi) ~ (0.6,20°). 
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Figure 13. Nondimensional time-averaged swimming speed and power consumption of 
squirmers with different initial position (/3/) and orientation (x/) with three-dimensional 

kinematics. 
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Figure 14. Two-dimensional trajectories, z{x), of pullers in a tube (left: a — 3; right: a = 5). 
Different combinations of the initial position (/3/) and pitching angle (^/) are chosen. The inset 
plots display the trajectories near the starting positions. 



trajectory, U, decrease clearly with xi but slowly with /3j. Larger values of xi and /3/ 
result in larger maximum off-axis distance, leading to higher hydrodynamic resistance 
from the boundaries and thus hindering locomotion. We also observe that U^ decreases 
with xi niore rapidly than U. As xi increases, the swimmer trajectory becomes more 
coiled, which significantly decreases the swimming velocity in the axial direction. We also 
note that the power consumption, V, increases with the the initial orientation, %/, but 
does not change significantly with /?/. 
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Figure 15. Equilibrium points, {rcYL,£,CYL}, of the puller in the tube as a function of the 
the dipole strength, a, and for a confinement a/R — 0.3. Blue dashed lines (a < 3.86) show 
the equilibrium point at (0, 0), corresponding to swimming in the centre of the tube and along 
its axis. For a > 3.86, the combination {rcYL,S,CYL} characterises the equilibrium state for 
swimming along a straight line with off-axis distance rcYL and orientation towards the wall 

^CYL- 



4.4. The trajectory of a puller inside the tube 
In this section, we study the trajectories of a puller swimmer (a > 0) in the tube. We 



first consider the case where the motion is restricted to the x — z plane, as in Sec. 4.2 In 
figure [14] we show the two-dimensional trajectories of pullers having dipole parameters 
of q; = 3 (left) and a = 5 (right), for different initial positions, j3j, and orientations, ^/. 
In both cases, the swimmers initially follow wavelike trajectories with decreasing magni- 
tude, and eventually settle along straight trajectories, displaying thus passive asymptotic 



stability (Yizhar & Richard 2009). The puller with a — 3 ends up swimming along the 



tube axis, with {rcYL,^CYL) = (0,0) as its equilibrium point (cylindrical coordinates 
are used here, and rcYL and S,cyl denote the off-axis distance and orientation of the 
cell respectively). In contrast, the puller with a — b swims parallel to the axis near the 
top or bottom wall depending on its initial position and orientation, thus its equilibrium 
point corresponds to swimming along an off-axis straight line. In that case, even though 
the trajectory is parallel to the tube axis, the swimmer remains slightly inclined towards 
the wall to offset the hydrodynamic repulsion from the wall. 

We further examine the coordinates of equilibrium points {rcYLiS^cYh) as a function 
of the dipole strength, a, in figure [15) For a below a critical value, ac ~ 3.86 for the 
confinement chosen here (a/R — 0.3), the equilibrium point is {rcYL,icYL) — (0,0) 
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Figure 16. Critical value of the dipole parameter, a^, for stable swimming of the puller in the 
tube centre as a function of the confinement, a/R. The approximate values of ac are denoted 
by the square symbols and the upper (resp. lower) limit of the error bar corresponds to the 
asymptotically-stable swimming motion away from the centre (resp. along the tube axis). 



denoted by the dashed blue line. For a > 3.86, the equilibrium point corresponds to 
swimming stably along a straight line with off-axis distance rcYL and orientation S,cyl , 
both of which grow with increasing a. The relationship between confinement, a/R, and 



the critical value ac is examined in figure 16 Determining precisely the value of ac 
is not possible due to the large computational cost so we report approximate values, 
with an upper (resp. lower) limit of the error bar corresponding to the asymptotically- 
stable swimming motion near the wall (resp. along the tube axis). The critical dipolar 
strength first increases with the confinement, reaching its maximum as a/R ~ 0.3, before 
decreasing. 

By starting with different combinations of a, j3i and xii ^^ obtain different three- 
dimensional trajectories for the puller. Some of these trajectories are illustrated in fig- 
ure 1171 Results similar to the two-dimensional simulations are obtained. For a below 
a critical value, pullers eventually swim along the tube axis indicating the equilibrium 
point {rcYLiS,cYL) = (OjO)- Foi' larger values of a, the equilibrium point corresponds 
to swimming motion with constant off-axis distance and orientation. Hydrodynamic in- 
teractions between the swimmer and the tube alone are responsible for such a passive 
stability, which could be of importance to guarantee, for example, robust steering of ar- 



tificial micro-swimmers in capillary tubes without on-board sensing and control (Yizhar 
fc Richard|[2009l ). 



We conclude this section by investigating in figure 18 the swimming speed of the puller 
along the stable trajectory and the dependency of its magnitude on the dipole strength, 
a. In the case of confinement a/R — 0.3, the swimming speed Ux is larger than that in 
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Figure 17. Three-dimensional trajectories of pullers in the tube with confinement a/R — 0.3. 
The red dashed line indicates the tube axis. Different combinations of the dipole strength, a, 
initial position, /?/, and initial orientation, x/i a-re chosen. Green trajectories correspond to 
pullers with one equilibrium point in the tube centre whereas blue ones are for pullers with 
equilibrium near the wall. 




Figure 18. Swimming speed, Ux (blue lower triangles), and power consumption, V (green upper 
triangles), of the puller as a function of the dipole parameter, a, both quantities being scaled by 
their corresponding values in free space (with a/R = 0.3). The critical value Oc = 3.86, shown 
by the cross, is the transition between stable swimming at the tube centre vs. stable swimming 
near the tube walls. The blue dashed line indicates the swimming speed in free space. 



free space as a is above a critical value (around 4 here) and it increases by about 16% 
as a = 5. This is an example of swimming microorganisms taking propulsive advantage 



from near-wall hydrodynamics, as discussed in previous analytical studies (Katz 1974 
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Figure 19. The orientation, e, of a puller swimming (red arrow) on the stable trajectory near 
the tube wall. Curved green arrows stand for the flow imposed at the surface of the swimmer, 
black (F^) and grey (Fp) dashed arrows for the hydrodynamic force while F_r is the repulsive 
force. 



Felderhof|2009 2010). In our case, as the squirmer is oriented into the wall, the direction 
of the wall-induced hydrodynamic force, Fn, resulting from flow being ejected on the 
side of the puller, is not normal to the wall but possesses a component in the swimming 
direction, as shown in figure |19| This force contributes thus to an additional propulsion 
and increases the swimming speed. 



4.5. The trajectory of a pusher inside the tube 

We next address the spherical pusher squirmer, with a negative force dipole, a. We find 
that the motion of the pushers inside the tube is unstable. The trajectories of pushers 



confined in the x — z plane (x/ = 0) are plotted in figure 20 for different combinations of 



dipole strength, initial position, and initial orientation. The pushers always execute wave- 
like motions with decreasing wavelengths and increasing amplitude, eventually crashing 
into the walls. Pushers and pullers display therefore very different swimming behaviours, 
a difference which stems from the opposite front-back asymmetry of the force dipole. 



4.6. Squirmers with normal surface velocity 

For the sake of completeness, we investigate in this section the dynamics of squirmers 
in the tube in the case where the squirming motion is induced by normal (instead of 
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Figure 20. Two-dimensional trajectories of pushers for different combinations of a, I3j, and 
^/. The black circles indicate the moment the swimmers make contact with the wall. 



tangential) surface velocity, modelled as 



UsAr(r) 



E 



^n{n + iy 



n=0 



(4.1) 



where An is the nth mode of the normal squirming velocity (Blake|[l971 ). In free space, 



the swimming velocity is Ugj^ = —Ai/3 (Blake|l971 ). For simplicity, we only consider the 



instantaneous kinematics of a squirmer with Ai = —1 and A^^i — 0, corresponding thus 
to Ugj^ ~ 1/3. The swimmer is located at (0, 0, — /? {R — a)), and is oriented in the positive 
X direction. We plot the axial velocity component, Ux, (scaled by Ugj^) together with the 
rotational velocity, Q.y, in figure 21 Both Ux and Vly are seen to increase monotonically 



with the confinement and eccentricity. This is in agreement with past mathematical 
analysis stating that microorganisms utilising transverse surface displacement speed up 
when swimming near walls ( Katz|[l974 ), between two walls ( Felderhof||2009 ) , or inside a 
tube ( |Felderhof|[2010l ). 

This increase (resp. decrease) of swimming speed in the tube of a squirmer with normal 
(resp. tangential) surface deformation can be related to the problem of micro-scale loco- 
motion in polymeric solutions. It is well known that actuated biological flagella generate 



drag-based thrust due to larger resistance to normal than to tangential motion (Lauga & 



Powers|[2009 1 . When swimming in polymer solutions, flagella undergoing motion normal 
to its shape push directly onto the neighbouring polymer network, whereas tangential 
motion barely perturb these micro obstacles ( Berg fc Turner|1979 Magariyama & Kudo 
2002 Nakamura et aL|2006 Leshansky|2009). In this case, the drag force increases more 
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Figure 21. Swimming velocity in the axial direction, Ux, and rotational velocity, Qy, of the 
squirmer with normal surface deformation with modes An = —5ni; Ux is scaled by the swimming 
speed in free space, Ugj^. The squirmer is located at (0, 0, — /? {R — a)) and oriented parallel to 
the axis. Different values of a/R are reported with maximum value oi P — 0.99. 



in the normal direction than in the tangential, resulting in larger swimming speeds (Berg 



fc Turner|[T979l [Magariyama fc Kudo||2002l [Nakamura eraI1[2006l |Leshansky|[2009l pu 



et al. 2011). Likewise, it was shown for a spherical squirmer that polymeric structures 



in the fluid always decrease the swimming speed in case of tangential surface deforma- 



tion (Leshansky 2009 Zhu et al. 2011 20121 but increase for normal deformation (Le- 



shansky]|2009 ) . The increase of swimming speed observed here in the case of a squirmer 
with normal surface deformation can similarly be attributed to the flow directly onto the 
tube wall. 



The value of rotational velocity, r2„, shown in figure 21 shows however that the squirmer 



rotates into the nearest wall, thus getting eventually trapped there. In order to avoid 
being trapped while at the same time taking advantage of the wall-induced enhanced 
propulsion, ideally swimmers should thus use a combination of tangential and normal 
deformation. 

Interestingly, a superposition of the neutral squirming mode (i3„ = (5„i, see q2]) with 
the first normal squirming mode (A„ = — ^ni) results in a special swimmer able to move 
without creating any disturbance in the surrounding fluid, characterised by a uniform 
squirming velocity of —1 everywhere on the body (in the co-moving frame), no body 
rotation, and a swimming speed equal to 1. This remains true regardless of the degree of 
conflnement as conflrmed by our numerical simulations. 
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Figure 22. Three-dimensional swimming of the neutral squirmer inside a torus-like curved 
tube. The dashed line indicates the circular axis of the torus with its radius, Rb = 20a, a is 
the azimuthal position of the squirmer, and R (cr) is the distance between the squirmer and the 
centre of the baseline circle. The wavelike motions for R (cr) — Rb and y{<j) are shown on the 
right. 

5. Swimming inside a curved tube 

In this final section, we investigate the squirmer motion inside a curved tube that is 
a part of a torus. The axis of the torus is a circle on the plane y = Q with its radius 
Rb = 20a. Trajectories of a neutral squirmer and a puller with the dipole strength a — 1 



are shown in Figs. 22 and 23 respectively. In both cases, the trajectory is displayed in 
both the X ~ z and y — z planes. The motion in the radial direction, represented by 
R {a) — Rb, is plotted as a function of the azimuthal position of the swimmer, tr, where 
R {a) is the distance between the cell and the centre of the circle. In both cases, the 
dynamics of swimmers initially starting aligned with the tube axis is wavelike. For the 
neutral squirmer, the wavelength and wave magnitude approach a constant value cr > tt 



(figure 22 right), indicating marginal stability of the motion. In contrast, for the puller, 



decaying waves are observed (figure 23 right), indicating passive asymptotic stability. 
As in the straight-tube case, pushers are unstable and crash into walls in finite time. 



6. Conclusion and outlook 

In this paper, a Boundary Element Method code was developed, validated, and used 
to present computations for the locomotion of model dilates inside straight and curved 
capillary tubes. We used the spherical squirmer as our model microorganism and studied 
the effect of confinement on the kinematics, energetics, and trajectories of the cell. We 
also investigated the stability of the swimming motion of squirmers with different gaits 
(neutral, pusher, puller). 
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Figure 23. Same as in figure [22] but for a puller witli a = 1. 



We found that tube confinement and near-wall swimming always decrease the swim- 
ming speed of a squirmer with tangential surface deformation for swimming parallel to 
the tube axis. In contrast, a swimmer with normal surface deformation improves its 
swimming speed by directly pushing against the surrounding tube wall. In both cases 
however, tube confinement and near-wall swimming always lead to additional viscous 
dissipation, thus increasing the power consumption. 

Focusing on swimming with tangential forcing, we then studied in detail the dynamics 
of neutral, puller, and pusher squirmers inside a straight tube. For a neutral squirmer, 
swimming motion on the tube axis is marginally stable and generically displays three- 
dimensional helical trajectories as previously observed experimentally for Paramecium 
cells. Importantly, these helical trajectories arise purely from hydrodynamic interactions 
with the boundaries of the tube. 

In the case of puller swimmers, their trajectories are wavelike with decreasing ampli- 
tude and increasing wavelength, eventually leading to stable swimming parallel with the 
tube axis with their bodies slightly oriented toward the nearest wall. The locations for 
these stable trajectories depend on the strength of the force dipole, a. Swimmers with 
weak dipoles (small a) swim in the centre of the tube while those with strong dipoles 
(large a) swim near the walls. The stable orientation of the swimmers makes an non-zero 
contribution of the wall-induced hydrodynamic forces in the direction of locomotion, thus 
leading to an increase of the swimming speed (although accompanied by an increase of 
the rate of viscous dissipation). In contrast, pushers are always unstable and crash into 
the walls of the tube in finite time. Similar results are observed for locomotion inside a 
curved tube. 

We envision that our study and general methodology could be useful in two specific 
cases. First, our results could help shed light on and guide the future design and ma- 
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neuverability of artificial small-scale swimmers inside small tubes and conduits. Second, 
the computational method could be extended to more complex, and biologically-relevant, 
geometries, to study for example the locomotion of flagellated bacteria or algae into con- 
fined geometries, as well as their hydrodynamic interactions with relevant background 
flows. It would be also interesting to relax some of our assumptions in future work, and 
address the role of swimmer geometry on their stability (we only considered the case of 
spherical swimmers in our paper) and quantify the role of noise and fluctuations on the 
asymptotic dynamics obtained here. 
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